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Implementing  Linear  Multistep  Formulas  for  Solving  DAEs 

INTRODUCTION 

Some  errors  were  discovered  in  the  numerical  experiments  of  the  report  "Implementing 
Linear  Multistep  Formulas  for  Solving  DAEs.”  (University  of  Illinois,  #UIUCDCS-R-85- 
1205,  April  1985.)  First,  due  to  a  coding  error,  the  error  estimate  involving  only  the 
differential  variable  was  incorrectly  implemented  for  the  full  index-2  problem. 

Also,  the  exact  value  of  the  derivative  of  the  constraint  function  was  incorrectly 
specified  for  the  full  index-2  problem  with  a  steep  gradient,  so  the  absolute  errors  in  u 
were  not  accurate.  Coincidentally,  the  incorrect  derivative  tended  to  mimic  the  behavior 
of  the  actual  derivative  over  the  interval  of  interest. 

The  following  is  a  complete  replacement  for  Section  14  of  the  report.  In  addition  to 
the  above  described  corrections,  it  includes  some  clarifications  and  corrections  to  the  text, 
many  more  tables  for  the  full  index-2  problem,  graphs  illustrating  the  difference  in 
efficiency  of  the  two  error  estimates  under  consideration,  and  simplified  results  for  the 


pendulum  problem. 


. 


14.  EXAMPLES 

We  now  consider  some  examples  from  the  literature. 


Problem  1 

z(t)  —  </>( t )  =  0 

z'(t )  —  u(t )  =  0 


Where  u0,  z0,  and  z0'  are  given,  and  <j)  is  a  given  function  defined  on  the  interval  of  interest. 
This  is  a  full  index  two  problem.  The  Jacobian  is  easily  computed: 
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So  X  =  (a4/A  -A3Af1A2)  1  =  a0 
The  error  estimate  is  given  by  Equation  (2.7) 
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And  we  have  e*  —  //2  and  e“  =  — rn  —  //4  +  o0//2/A.  Hence  2  is  computed  as  accurately  as  is 
computed  and  errors  in  u  depend  on  the  local  truncation  error  in  approximating  z1  as  well  as  on 
how  accurately  the  equations  2—  (j)[t)  =  0  are  solved. 


We  applied  the  double— precision  version  of  the  code  DASSL  of  Petzold  (1982b)  to  this  problem. 
We  were  most  interested  in  the  behavior  of  the  code  when  the  function  </>  has  a  steep  gradient  at 
some  point  in  the  interval  of  integration,  since  the  theory  predicts  that  standard  ODE  error  esti¬ 
mates  will  be  especially  unreliable  in  that  case.  As  a  control,  though,  we  also  examined  the 
integrator’s  performance  when  solving  a  problem  with  a  relatively  smooth  driving  function  (j). 

The  smooth  function  was 

4>(t)  =  5tn(27T<). 


And  the  problem  was  solved  from  t  =0.0  to  t  =1.0,  with  exact  initial  values  u(0)=27T  and 

2(0)=0.0. 

Although  DASSL  accepts  arbitrary  DAEs  of  the  form  F(y,y',t)  =  0,  it  is  designed  to  handle 
problems  of  index  no  greater  than  one.  We  modified  DASSL  to  control  only  errors  in  the 
differential  variables  {z)  by  assuming  that  the  equations  are  so  arranged  that  the  differential  vari¬ 
ables  are  the  first  components  of  the  solution  and  the  state  variables  (u)  the  last.  The  error  con¬ 
trol  on  truncation  error  now  operates  on  just  the  first  “NDIFF”  components  of  y. 

Table  14.1a  gives  the  exact  errors  in  the  computed  solution  when  NDIFF=2.  This  means 
the  truncation  errors  in  both  variables  z  and  u  were  controlled.  For  various  tolerances  we  give 
the  exact  errors  e(z)  and  e(u)  at  the  endpoint.  For  this  problem,  the  errors  were  of  similar  mag¬ 
nitude  throughout  the  interval  of  integration. 

Because  DASSL  does  not  re-evaluate  the  jacobian  at  each  step,  even  the  algebraic  equation 
*(0  =  is  not  solved  exactly,  so  there  is  some  error  in  z.  Furthermore,  DASSL  computes 
values  at  the  endpoint  by  interpolation,  so  additional  interpolation  errors  are  introduced  at  t=l. 
Interpolation  errors  in  u  are  not  serious  because  that  variable  already  has  relatively  large 
integration  errors  at  the  endpoint,  but  z  is  computed  quite  accurately  throughout  the  interval,  so 
the  error  in  z  at  t=l  will  reflect  the  size  of  the  step,  and  therefore  the  tolerance. 


Table  14.1a.  Exact  errors  at  f=l,  NDIFF=2. 


TOL 

efz)  fat  t=1.0) 

o 

r-H 

II 

o 3 

^CD 

0.10000e+00 

0.20416e-02 

0.14571e+00 

0.50000e-01 

0.17564e-03 

0.38709e-01 

0.25000e-01 

0.30873e-04 

0.29924e-01 

0.12500e-01 

0.10716e-04 

0.21103e-01 

0.62500e-02 

0.14568e-04 

0.42532e-02 

0.31250e-02 

0.44214e-04 

0.50268e-02 

0.15625e-02 

0.99309e-07 

0.43658e-02 

0.78125e-03 

0.22125e-05 

0.52767e-02 

0.39063e-03 

0.13470e-05 

0.14117e-02 

When  NDIFF=1,  only  the  errors  in  variable  z  were  controlled.  The  errors  //j  and  //2  are 
controlled  because  DASSL  controls  the  iteration  error  in  the  corrector  solution. 


Table  14.1b.  Exact  errors  at  t=  1,  NDIFF=1. 


TOL 

efz)  (at  t=  1.0) 

e(u)  (at  t=1.0) 

0.10000e+00 

0.40835e-01 

0.88190e+00 

0.50000e-01 

0.21189e-02 

0.22352e+00 

0.25000e-01 

0.10708e-02 

0.66365e-01 

0.12500e-01 

0.34440e-03 

0.47012e-01 

0.62500e-02 

0.31464e-03 

0.28675e-01 

0.31250e-02 

0.27920e-04 

0.21890e-01 

0.15625e-02 

0.86727e-05 

0.39854e-02 

0.78125e-03 

0.48029e-04 

0.61423e-03 

0.39063e-03 

0.22780e-04 

0.27375e-03 

The  work  performed  at  a  given  tolerance  depends  on  whether  or  not  error  in  both  variables 
is  being  controlled.  Table  14.1c  summarizes  the  corresponding  amount  of  work  performed  in 
each  case. 


Table  14.1c.  Work  performed  for  various  tolerances. 


TOL 

NDIFF 

No.  Steps 

0.10000e+00 

2 

23 

1 

14 

0.50000e-01 

2 

32 

1 

16 

0.25000e-01 

2 

70 

1 

18 

0. 12500e-01 

2 

76 

1 

19 

0.62500e-02 

2 

130 

1 

22 

0.31250e-02 

2 

388 

1 

24 

0. 15625e-02 

2 

109 

1 

27 

0.78125e-03 

2 

1483 

1 

30 

0.39063e-03 

2 

90 

1 

34 

Thus  while  there  may  be  some  decrease  in  accuracy  when  the  error  control  is  removed  from 
u,  there  is  a  tremendous  decrease  in  the  amount  of  work  performed. 

In  figure  14.1,  we’ve  plotted  the  actual  error  in  u  at  t=  1  versus  the  number  of  steps  (these 
results  were  generated  by  running  the  code  for  various  tolerances).  A  logarithmic  scale  is  used 
because  of  the  great  variations  in  magnitude  involved.  The  graph  clearly  demonstrates  the  great 
superiority  of  the  selective  error  control — even  for  this  simple  problem  with  a  smooth  gradient. 


Figure  14.1.  Accuracy  in  u  vs.  Work  Performed,  smooth  constraint. 
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Number  of  Succesful  Steps 


The  second  function  is  given  by 
d>(t)  =  sin(2iTt)-\-tanh(30t— 15). 

This  function  has  a  steep  gradient  in  the  middle  of  the  interval  (0.0,  1.0).  Again  the  exact  initial 
values  given  below  were  supplied  to  the  code. 

u(0)=2/T+30sec/i(— 15)  and  z(0)=tanh(— 15). 

As  has  been  shown,  a  BDF  code  using  an  error  estimate  based  on  the  difference  between  the  pred¬ 
ictor  and  the  corrector  cannot  be  expected  to  compute  accurate  approximations  to  the  algebraic 
variable  in  the  neighborhood  of  a  steep  gradient  (see  Petzold  (1982a)). 

Table  14.2a  gives  the  actual  absolute  and  relative  error  in  u  and  absolute  error  in  z  at 
selected  points  in  the  interval  for  a  tolerance  of  0.1.  As  before,  interpolation  errors  exaggerate 
inaccuracy  in  z  at  the  endpoint.  It  is  interesting  to  note  the  points  at  which  DASSL  did  not  re¬ 
evaluate  the  jacobian:  these  appear  as  a  large  error  in  z  relative  to  machine  round-off. 


Table  14.2a.  Errors  at  selected  points,  NDIFF=2 


TOL 

Time 

efu) 

rel.  e(u) 

e(z) 

0.10000e+00 

0.10000e-02 

0.82690e-04 

0.13161e-04 

0.71999e-ll 

0.31000e-01 

0.50699e-01 

0.81574e-02 

0.21898e-09 

0.63000e-01 

0.20376e+00 

0.33953e-01 

0.41852e-09 

0.42567e+00 

0.48183e+00 

0.10315e+00 

0.27367e-08 

0.43984e+00 

0.12149e+01 

0.33534e+00 

0.29434e-08 

0.45260e+00 

0.33019e+01 

0.22986e+01 

0.31152e-08 

0.46779e+00 

0.11247e+02 

0.17722e+01 

0.32993e-08 

0.49417e+00 

0.94514e+03 

0.46146e+02 

0.35557e-08 

0.51 178e+00 

0.20219e+03 

0.91239e+01 

0.36747e-08 

0.52749e+00 

0.16478e+02 

0.12992e+01 

0.37414e-08 

0.54583e+00 

0.12234e+01 

0.81111e+00 

0.37682e-08 

0.56928e+00 

0.11802e+00 

0.32419e-01 

0.37166e-08 

0.58727e+00 

0.91314e-01 

0.19742e-01 

0.36085e-08 

0.94437e+00 

0.23498e+00 

0.38282e-01 

0.63883e-08 

0.10000e+01 

0.12049e+00 

0.19552e-01 

0.46832e-03 

Keeping  NDIFF=2,  we  see  that  at  smaller  tolerances,  the  errors  in  z  quickly  become  negligi¬ 
ble  throughout  the  interval.  Tables  14.2b— i  give  exact  errors  in  u  and  z  at  selected  tolerances 
and  for  selected  intermediate  steps. 


Table  14.2b.  Errors  at  selected  points,  NDIFF=2 


TOL 

Time 

e(u) 

rel.  efu) 

efz) 

0.50000e-01 

0. 10000e-02 

Q.82690e-04 

0.13161e-04 

0.71999e-ll 

0.31000e-01 

0.50699e-01 

0.81574e-02 

0.21898e-09 

0.63000e-01 

0.20376e+00 

0.33953e-01 

0.41852e-09 

0.4 1601e+00 

0.40385e+0Q 

0.80054e-01 

0.25878e-08 

0.43376e+00 

0.52946e+0Q 

0.13457e+00 

0.28564e-08 

0.44973e+00 

0.19038e+01 

0.13824e+01 

0.30779e-08 

0.45623e+00 

0.27674e+01 

0.21707e+01 

0.31614e-08 

0.47639e+00 

0.33544e+02 

0.30555e+01 

0.33924e-08 

0.49839e+00 

0.12744e+05 

0.54977e+03 

0.35882e-08 

0.51939e+00 

0.56406e+02 

0.33834e+01 

0.371 19e-08 

0.53708e+00 

0.58891e+01 

0.13672eT01 

0.37625e-08 

0.56084e+00 

0.47232e+00 

0.15662e+00 

0.37466e-08 

0.58573e+00 

0.10930e+Q0 

0.23888e-01 

0.36201e-08 

0.93654e+00 

0.16465e+00 

0.29268e-01 

0.62140e-08 

0.95613e+00 

0.12042e+00 

0.20322e-01 

0.66243e-08 

0.99533e+00 

0.85874e-01 

0.13863e-01 

0.71633e-08 

0.10000e+01 

0.83201e-01 

0.13420e-01 

0.19849e-03 

Table  14.2c.  Errors  at  selected  points,  NDIFF=2 


TOL 

Time 

- .  .?(u.)  . 

rel.  efu) 

efz) 

0.25000e-01 

0.10000e-02 

0.82690e-04 

0.13161e-04 

0.71999e-ll 

0.31000e-01 

0.50699e-01 

0.81574e-02 

0.21898e-09 

0.63000e-01 

Q.20376e+00 

0.33953e-01 

0.41852e-09 

0.41395e+00 

0.85246e-01 

0.17843e-01 

0.25553e-08 

0.42834e+00 

0.31094e+00 

0.72420e-01 

0.27767e-08 

0.44871e+00 

0.11725e+01 

0.1 1098e+01 

0.30644e-08 

0.46171e+00 

0.51744e+01 

0.14241e+01 

0.32286e-08 

0.48160e+00 

0.67318e+02 

0.43581e+01 

0.34443e-08 

0.5051  le+00 

0.12354e+04 

0.52046e+02 

0.36349e-08 

0.52293e+00 

0.34864e+02 

0.26404e+01 

0.37261e-08 

0.54221e+00 

0.28519e+01 

0.12197e+01 

0.37674e-08 

0.56163e+00 

0.15783eT00 

0.55378e-01 

0.37443e-08 

0.58374e+00 

0.76270e-01 

0.16735e-01 

0.36344e-08 

0.94597e+00 

0.17085e+00 

0.29695e-01 

0.64222e-08 

0.96888e+00 

0.10901e+00 

0.18005e-01 

0.68431e-08 

0.99180e+00 

0.44914e  -01 

0.72094e-02 

0.71315e-08 

0.10000e+01 

0.16156e-01 

0.25779e-02 

0.15494e-03 

Table  14. 2d.  Errors  at  selected  points,  NDIFF=2 


TOL 

Time 

eM 

rel.  efu) 

e(z) 

0.12500e-01 

0.10000e-02 

0.82690e-04 

0.13161e-04 

0.71999e-ll 

0.31000e-01 

0.50699e-01 

0.81574e-02 

0.21898e-09 

0.63000e-01 

0.20376e+00 

0.33953e-01 

0.41852e-09 

0.41938e+00 

0.11148e+00 

0.24030e-01 

0.26403e-08 

0.43774e+00 

0.45095e+00 

0.13858e+00 

0.29137e-08 

0.45119e+00 

0.14612e+01 

0.507 10e+01 

0.30970e-08 

0.47217e+00 

0.18815e+02 

0.20258e+01 

0.33478e-08 

0.4927  le+00 

0.58896e+03 

0.26877e+02 

0.35438e-08 

0.51244e+00 

0.17948e+03 

0.89410e+01 

0.36783e-08 

0.53110e+00 

0.11671e+02 

0.14325e+01 

0.37512e-08 

0.55096e+00 

0.10138e+01 

0.13228e+01 

0.37653e-08 

0.57 142e+00 

0.94616e-02 

0.23822e-02 

0.37069e-08 

0.59320e+00 

0.76077e-02 

0.15922e-02 

0.35594e-08 

0.95315e+00 

0.92776e-01 

0.15671e-01 

0.65675e-08 

0.98258e+00 

0.17608e-01 

0.281 13e-02 

0.70322e-08 

0.99669e+00 

0.14015e-01 

0.22360e-02 

0.71746e-08 

0.10000e+01 

0.12085e-01 

0.19270e-02 

0.77253e-05 

Table  14. 2e.  Errors  at  selected  points,  NDIFF=2 


TOL 

Time 

efu) 

rel.  efu) 

efz) 

0.62500e-02 

0.10000e-02 

0.82690e-04 

0.13161e-04 

0.71999e-ll 

0.31000e-01 

0.50699e-01 

0.81574e-02 

0.21898e-09 

0.47000e-01 

0.81720e-01 

0.13412e-01 

0.32375e-09 

0.61400e-01 

0.98946e-01 

0.16713e-01 

0.40959e-09 

0.40322e+00 

0.26159e-02 

0.54598e-03 

0.23828e-08 

0.42064e+00 

0.14144e+00 

0.30644e-01 

0.26598e-08 

0.44235e+00 

0.69610e+00 

0.27338e+00 

0.29783e-08 

0.46006e+00 

0.41717e+01 

0.14047e+01 

0.32087e-08 

0.48068e+00 

0.58521e+02 

0.38471e+01 

0.34355e-08 

0.50056e+00 

0.10475e+06 

0.43914e+04 

0.36040e-08 

0.5202  le+00 

0.51070e+02 

0.34027e+01 

0.37154e-08 

0.54111e+00 

0.34945e+01 

0.13493e+01 

0.37667e-08 

0.56137e+00 

0.297 10e+00 

0.10099e+00 

0.37451e-08 

0.58160e+00 

0.62287e-02 

0.13665e-02 

0.36491e-08 

0.96259e+00 

0.10356e-01 

0.16919e-02 

0.67400e-08 

0.10000e+01 

0.58995e-02 

0.93806e-03 

0.20047e-05 

Table  14. 2f.  Errors  at  selected  points,  NDIFF=2 


TOL 

Time 

efu) 

rel.  efu) 

efz) 

0.31250e-02 

0.70337e-03 

0.409 13e-04 

0.65116e-05 

0.50642e-ll 

0.21805e-01 

0.25141e-01 

0.40229e-02 

0.15552e-09 

0.45283e-01 

0.55608e-01 

0.91368e-02 

0.31293e-09 

0.40034e+00 

0.22314e-01 

0.46406e-02 

0.23356e-08 

0.42182e+00 

0.69883e-01 

0.15571e-01 

0.26780e-08 

0.44099e+00 

0.40439e+00 

0.15734e+00 

0.29594e-08 

0.46022e+00 

0.41293e+01 

0.13055e+01 

0.32107e-08 

0.48168e+00 

0.67342e+02 

0.41403e+01 

0.34451e-08 

0.50029e+00 

0.40551e+06 

0.17116e+05 

0.36021e-08 

0.52355e+00 

0.32129e+02 

0.25334e+01 

0.37284e-08 

0.54046e+00 

0.37942e+01 

0.13268e+01 

0.37663e-08 

0.56190e+00 

0.28209e+00 

0.93327e-01 

0.37435e-08 

0.58052e+00 

0.28373e-01 

0.62334e-02 

0.36561e-08 

0.94508e+00 

0.41294e-01 

0.70331e-02 

0.64034e-08 

0.96563e+00 

0.33493e-01 

0.54872e-02 

0.67910e-08 

0.98042e+00 

0.20156e-01 

0.32429e-02 

0.70057e-08 

0.10000e+01 

0.50726e-02 

0.80667e-03 

0.11121e-05 

Table  14. 2g.  Errors  at  selected  points,  NDIFF=2 


TOL 

Time 

efu) 

rel.  efu) 

efz) 

0.15625e-02 

0.35169e-03 

0.10234e-04 

0.16288e-05 

0.25321e-l  1 

0.22156e-01 

0.25628e-01 

0.41018e-02 

0.15798e-09 

0.45417e-01 

0.33109e-01 

0.546 16e -02 

0.31378e-09 

0.64127e-01 

0.29942e-02 

0.51776e-03 

0.42474e-09 

0.40139e+00 

0.12017e-01 

0.25024e-02 

0.23528e-08 

0.42291e+00 

0.98865e-01 

0.22166e-01 

0.26946e-08 

0.44019e+00 

0.39686e+00 

0.14512e+00 

0.29481e-08 

0.46009e+00 

0.40830e+01 

0.13235e+01 

0.32091e-08 

0.48010e+00 

0.53540e+02 

0.35477e+01 

0.34298e-08 

0.50137e+00 

0.17628e+05 

0.74378e+03 

0.36098e-08 

0.52003e+00 

0.52396e+02 

0.34592e-f01 

0.37146e-08 

0.54005e+00 

0.39737e+01 

0.12967e+01 

0.37659e-08 

0.56022e+00 

0.33091e+00 

0.12028e+00 

0.37483e-08 

0.58032e+00 

0.13101e-01 

0.28929e-02 

0.36574e-08 

0.92892e+00 

0.19158e-02 

0.33796e-03 

0.60321e-08 

0.97005e+00 

0.3601 le-02 

0.58376e-03 

0.68611e-08 

0.10000e+01 

0.66383e-02 

0.10576e-02 

0.18954e-04 

Table  14. 2h.  Errors  at  selected  points,  NDIFF=2 


TOL 

Time 

e(u)  . 

rel.  e(u] 

efz) 

0.78125e-03 

0.17584e-03 

0.25638e-05 

0.40805e-06 

0.12661e-ll 

0.21769e-01 

0.12579e-01 

0.20168e-02 

0.15528e-09 

0.41439e-01 

0.73309e-02 

0.12060e-02 

0.28830e-09 

0.60358e-01 

0.21893e-01 

0.37369e-02 

0.40370e-09 

0.40149e+00 

0.16419e-03 

0.34273e-04 

0.23544e-08 

0.42022e+00 

0.55828e-01 

0.12274e-01 

0.26533e-08 

0.44000e+00 

0.36839e+00 

0.13424e+00 

0.29456e-08 

0.46002e+00 

0.40354e+01 

0.13174e+01 

0.32082e-08 

0.48007e+00 

0.53204e+02 

0.35194e+01 

0.34294e-08 

0.50006e+00 

0.87531e+07 

0.36908e+06 

0.36004e-08 

0.52106e+00 

0.45258e+02 

0.31458e+01 

0.37189e-08 

0.54152e+00 

0.33357e+01 

0.13879e+01 

0.37670e-08 

0.56001e+00 

0.35730e+00 

0.13060e+00 

0.37488e-08 

0.58051e+00 

0.18964e-01 

0.41753e-02 

0.36562e-08 

0.95760e+00 

0.64752e-02 

0.10671e-02 

0.66515e-08 

0.98590e+00 

0.24461e-02 

0.39068e-03 

0.70707e-08 

0.10000e+01 

0.89851e-03 

0.14298e-03 

0.50553e-05 

Table  14. 2i.  Errors  at  selected  points,  NDIFF=2 


TOL 

Time 

efu) 

rel.  efu) 

efz) 

0.39063e-03 

0.87922e-04 

0.64636e-06 

0.10287e-06 

0.63302e-12 

0.22420e-01 

0.26007e-02 

0.41788e-03 

0.15982e-09 

0.40971e-01 

0.44280e-02 

0.72822e-03 

0.28527e-09 

0.62864e-01 

0.41193e-02 

0.70980e-03 

0.41777e-09 

0.40389e+00 

0.91883e-02 

0.19130e-02 

0.23937e-08 

0.42010e+00 

0.39300e-01 

0.86620e-02 

0.26515e-08 

0.44005e+00 

0.36666e+00 

0.13412e+00 

0.29462e-08 

0.46002e+00 

0.40300e+01 

0.13136e+01 

0.32082e-08 

0.48002e+00 

0.52861e+02 

0.34984e+01 

0.34290e-08 

0.50025e+00 

0.53581e+06 

0.22593e+05 

0.36018e-08 

0.52038e+00 

0.49882e+02 

0.33576e+01 

0.37161e-08 

0.54102e+00 

0.35462e+01 

0.13550e+01 

0.37667e-08 

0.56113e+00 

0.31320e+00 

0.10758e+00 

0.37458e-08 

0.58111e+00 

0.26838e-01 

0.58687e-02 

0.36523e-08 

0.94706e+00 

0.36128e-03 

0.60838e-04 

0.64451e-08 

0.96951e+00 

0.427 16e-03 

0.69256e-04 

0.68528e-08 

0.98972e+00 

0.22773e-03 

0.36321e-04 

0.71111e-08 

0.10000e+01 

0.15385e-03 

0.24487e-04 

0.10039e-05 

We  now  consider  the  results  for  the  case  when  only  the  error  in  variable  z  is  controlled  (i.e. 
NDIFF=1).  Table  14.3a  gives  the  results  for  TOL=0.1. 


Table  14.3a.  Errors  at  selected  points,  NDIFF=1 


TOL 

Time 

. — ... 

rel.  efu) 

efz) 

0.10000e+00 

0.10000e-02 

0.82690e-04 

0.13161e-04 

0.71999e-ll 

0.31000e-01 

0.11405e-03 

0.18502e-04 

0.21898e-09 

0.63000e-01 

0.18741e-02 

0.32338e-03 

0.41852e-09 

0.44543e+00 

0.34267e+01 

0.77194e+00 

0.30204e-08 

0.48269e+00 

0.84354e+02 

0.75546e+01 

0.34548e-08 

0.51995e+00 

0.39766e+02 

0.14002e+01 

0.37143e-08 

0.53873e+00 

0.70978e+01 

0.57242e+01 

0.37646e-08 

0.55751e+00 

0.23621e+01 

0.56534e+00 

0.37549e-08 

0.58641e+00 

0.79843e+00 

0.20465e+00 

0.36150e-08 

0.92342e+00 

0.50469e+00 

0.99642e-01 

0.58936e-08 

0.10000e+01 

0.36701e+00 

0.55188e-01 

0.63112e-02 

At  smaller  tolerances,  we  see  the  same  expected  effects  as  before,  tables  14.3b-i  give  results 
for  NDIFF=1. 


Table  14.3b.  Errors  at  selected  points,  NDIFF=1 


TOL 

Time 

eM 

rel.  e(u) 

e(z) 

0.50000e-01 

0.10000e-02 

0.82690e-04 

0.13161e-04 

0.71999e-ll 

0.31000e-01 

0.11405e-03 

0.18502e-04 

0.21898e-09 

0.63000e-01 

0.18741e-02 

0.32338e-03 

0.41852e-09 

0.43270e+00 

0.17078e+01 

0.32561e+00 

0.28410e-08 

0.46455e+00 

0.10349e+02 

0.52750e+01 

0.32622e-08 

0.48561e+00 

0.12739e+03 

0.71845e+01 

0.34821e-08 

0.50667e+00 

0.70481e+03 

0.25451e+02 

0.36448e-08 

0.51985e+00 

0.49755e+02 

0.25833e+01 

0.37139e-08 

0.53303e+00 

0.53046e+01 

0.49153e+00 

0.37555e-08 

0.55398e+00 

0.70844e+00 

0.59640e+01 

0.37614e-08 

0.57492e+00 

0.34625e+00 

0.89152e-01 

0.36892e-08 

0.59587e+00 

0.14905e+00 

0.32090e-01 

0.35352e-08 

0.91005e+00 

0.32679e+00 

0.58014e-01 

0.55334e-08 

0.99383e+00 

0.51684e+00 

0.76058e-01 

0.71502e-08 

0.10000e+01 

0.52498e+00 

0.77111e-01 

0.17552e-02 

Table  14.3c.  Errors  at  selected  points,  NDIFF=1 


TOL 

Time 

e(u) 

rel.  e(u) 

e(z) 

0.25000e-01 

0.10000e-02 

0.82690e-04 

0.13161e-04 

0.71999e-ll 

0.31000e-01 

0.11405e-03 

0.18502e-04 

0.21898e-09 

0.63000e-01 

0.18741e-02 

0.32338e-03 

0.41852e-09 

0.44700e+00 

0.29260e+01 

0.85275e+00 

0.30416e-08 

0.45621e+00 

0.37404e+00 

0.10236e+00 

0.31612e-08 

0.47463e+00 

0.25064e+02 

0.21655e+01 

0.33740e-08 

0.49305e+00 

0.64875e+03 

0.26423e+02 

0.35466e-08 

0.50727e+00 

0.59313e+03 

0.27234e+02 

0.36485e-08 

0.5357  le+00 

0.74251e+01 

0.16406e+01 

0.37604e-08 

0.54993e+00 

0.17632e+01 

0.16148e+01 

0.37662e-08 

0.57553e+00 

0.16309eT00 

0.39721e-01 

0.36859e-08 

0.58833e+00 

0.12958e+00 

0.28138e-01 

0.36002e-08 

0.90829e+00 

0.42391e-01 

0.7981 7e-02 

0.54837e-08 

0.95949e+00 

0.30874e-01 

0.50517e-02 

0.66857e-08 

O.lOOOOe  mi 

0.37829e-02 

0.60170e-03 

0.46397e-03 

Table  14.3d.  Errors  at  selected  points,  NDIFF=1 


TOL 

Time 

e( u) 

rel.  efu) 

e(z) 

0.12500e-01 

0.10000e-02 

0.82690e-04 

0.13161e-04 

0.71999e-ll 

0.31000e-01 

0.11405e-03 

0.18502e-04 

0.21898e-09 

0.63000e-01 

0.18741e-02 

0.32338e-03 

0.41852e-09 

0.44700e+00 

0.29260e+01 

0.85275e+00 

0.30416e-08 

0.44897e+00 

0.13114e+01 

0.12009e+01 

0.30679e-08 

0.46791e+00 

0.12593e+02 

0.24097e+01 

0.33007e-08 

0.48988e+00 

0.28967e+03 

0.14565e+02 

0.35200e-08 

0.51318e+00 

0.15437e+03 

0.71130e+01 

0.36822e-08 

0.53065e+00 

0.11202e+02 

0.11677e+01 

0.37501e-08 

0.55037e+00 

0.16973e+01 

0.14031e+01 

0.37659e-08 

0.57008e+00 

0.38652e+00 

0.91493e-01 

0.37131e-08 

0.58979e+00 

0.40656e-01 

0.86229e-02 

0.35884e-08 

0.92640e+00 

0.14668e+00 

0.25421e-01 

0.59695e-08 

0.97719e+00 

0.19653e+00 

0.30634e-01 

0.69637e-08 

0.10000e+01 

0.20624e+00 

0.31781e-01 

0.19205e-02 

Table  14. 3e.  Errors  at  selected  points,  NDIFF=1 


TOL 

Time 

efu) 

rel.  efu) 

e(z) 

0.62500e-02 

0.10000e-02 

0.82690e-04 

0.13161e-04 

0.71999e-ll 

0.31000e-01 

0.11405e-03 

0.18502e-04 

0.21898e-09 

0.63000e-01 

0.18741e-02 

0.32338e-03 

0.41852e-09 

0.41628e+00 

0.41275e+00 

0.81797e-01 

0.25920e-08 

0.44232e+00 

0.95904e+00 

0.34032e+00 

0.29778e-08 

0.46066e+00 

0.46999e+01 

0.15592e+01 

0.32159e-08 

0.47900e+00 

0.45690e+02 

0.31750e+01 

0.34188e-08 

0.49501e+00 

0.12966e+04 

0.51558e+02 

0.35623e-08 

0.50867e+00 

0.40497e+03 

0.17945e+02 

0.36570e-08 

0.52320e+00 

0.32522e+02 

0.23065e+01 

0.37271e-08 

0.54729e+00 

0.19349e+01 

0.18318e+03 

0.37678e-08 

0.56794e+00 

0.12777e+00 

0.36672e-01 

0.37222e-08 

0.59547e+00 

0.12894e+00 

0.27644e-01 

0.35389e-08 

0.92585e+00 

0.44557e-01 

0.78750e-02 

0.59555e-08 

0.98091e+00 

0.327 13e-01 

0.52168e-02 

0.70118e-08 

0.10000e+01 

0.27063e-01 

0.42888e-02 

0.98061e-04 

Table  14. 3f.  Errors  at  selected  points,  NDIFF=1 


TOL 

Time 

e(u) 

rel.  efu) 

.  . e(z)  . 

0.31250e-02 

0.70337e-03 

0.40913e-04 

0.651 16e-05 

0.50642e-ll 

0.21805e-01 

0.27939e-04 

0.44886e-05 

0.15552e-09 

0.44312e-01 

0.4607 le-03 

0.76268e-04 

0.30676e-09 

0.89328e-01 

0.73933e-02 

0.13918e-02 

0.54449e-09 

0.40444e+00 

0.12517e+00 

0.25448e-01 

0.24026e-08 

0.4331  le+00 

0.41261e+00 

0.10582e+00 

0.28471e-08 

0.43833e+00 

0.15389e+00 

0.60486e-01 

0.29221e-08 

0.45921e+00 

0.36693e+01 

0.13606e+01 

0.31983e-08 

0.48009e+00 

0.53128e+02 

0.34480e+01 

0.34296e-08 

0.49719e+00 

0.41968e+04 

0.17843e+03 

0.35792e-08 

0.51719e+00 

0.80268e+02 

0.47771e+01 

0.37021e-08 

0.53931e+00 

0.41760e+01 

0.11694e+01 

0.37652e-08 

0.55689e+00 

0.45198e-f00 

0.21391e+00 

0.37562e~08 

0.57447e+00 

0.74719e-01 

0.17481e-01 

0.36916e-08 

0.90850e+00 

0.96116e-02 

0.18261e-02 

0.54896e-08 

0.97179e+00 

0.19992e-01 

0.32429e-02 

0.68872e-08 

0.10000e+01 

0.20518e-01 

0.32763e-02 

0.18892e-03 

Table  14. 3g.  Errors  at  selected  points,  NDIFF=1 


TOL 

Time 

e(u) 

rel.  e(u) 

e(z) 

0.15625e-02 

0.35169e-03 

0.10234e-04 

0.16288e-05 

0.25321e-ll 

0.22156e-01 

0.28881e-04 

0.46414e-05 

0.15798e-09 

0.44664e-01 

0.46819e-03 

0.77554e-04 

0.30900e-09 

0.89680e-01 

0.7451  le-02 

0.14047e-02 

0.54587e-09 

0.40205e+00 

0.71414e-01 

0.14684e-01 

0.23636e-08 

0.42193e+00 

0.14782e+00 

0.32413e-01 

0.26797e-08 

0.43660e+00 

0.34677e+00 

0.10377e+00 

0.28973e-08 

0.45126e+00 

0.16624e+01 

0.36326e+01 

0.30979e-08 

0.47074e+00 

0.15250e+02 

0.17163e+01 

0.33322e-08 

0.48708e+00 

0. 16389e+03 

0.82830e+01 

0.34954e-08 

0.50495e+00 

0.13213e+04 

0.57267e+02 

0.36339e-08 

0.52632e+00 

0.22365e+02 

0.20895e+01 

0.37378e-08 

0.54234e+00 

0.29989e+01 

0.14367e+01 

0.37675e-08 

0.56905e+00 

0.1 1877e+00 

0.30829e-01 

0.37176e-08 

0.5904  le+00 

0.40426e-01 

0.84161e-02 

0.35833e-08 

0.92159e+00 

0.23197e-01 

0.42079e-02 

0.58463e-08 

0.99238e+00 

0.32575e-01 

0.52175e-02 

0.7 1369e-08 

0.  10000e+01 

0.32216e-01 

0.51538e-02 

0.82566e-04 

Table  14. 3h.  Errors  at  selected  points,  NDIFF=1 


TOL 

Time 

e(u) 

rel.  e(u) 

e(z) 

0.78125e-03 

0.17584e-03 

0.25638e-05 

0.40805e-06 

0.12661e-ll 

0.22332e-01 

0.29351e-04 

0.47178e-05 

0.15921e-09 

0.44840e-01 

0.47193e-03 

0.78198e-04 

0.31012e-09 

0.89856e-01 

0.74801e-02 

0.14112e-02 

0.54656e-09 

0.41344e+00 

0.60675e-01 

0.12737e-01 

0.25473e-08 

0.42709e+00 

0.13108e+00 

0.31086e-01 

0.27580e-08 

0.44073e+00 

0.54101e+00 

0.19582e+00 

0.29558e-08 

0.45848e+00 

0.34125e+01 

0.14641e+01 

0.31894e-08 

0.48333e+00 

0.86598e+02 

0.49315e+01 

0.34609e-08 

0.49837e+00 

0.12487e+05 

0.52701e+03 

0.35881e-08 

0.51704e+00 

0.82034e+02 

0.48192e-f01 

0.37014e-08 

0.53968e+00 

0.41118e+01 

0.12550e+01 

0.37656e-08 

0.55963e+00 

0.36862e+00 

0.13817e+00 

0.37498e-08 

0.58623e+00 

0.75778e-01 

0.15880e-01 

0.36163e-08 

0.59813e+00 

0.43051e-01 

0.89046e-02 

0.35135e-08 

0.93017e+00 

0.41206e-02 

0.72496e-03 

0.60628e-08 

0.97301e+00 

0.53161e-02 

0.85913e-03 

0.69052e-08 

0.10000e+01 

0.57782e-02 

0.92048e-03 

0.39493e-04 

Table  14.31.  Errors  at  selected  points,  NDIFF=1 


TOL 

Time 

e(u) 

rel.  e(u) 

efz) 

0.39063e-03 

0.87922e-04 

0.64636e-06 

0.10287e-06 

0.63302e-12 

0.22420e-01 

0.29587e-04 

0.47560e-05 

0.15982e-09 

0.44928e-01 

0.47380e-03 

0.78521e-04 

0.31068e-09 

0.89944e-01 

0.74945e-02 

0.14144e-02 

0.54691e-09 

0.40806e+00 

0.31157e-01 

0.64797e-02 

0.24613e-08 

0.42330e+00 

0.70094e-01 

0.15896e-01 

0.27006e-08 

0.43570e+00 

0.22904e+00 

0.68098e-01 

0.28844e-08 

0.45925e+00 

0.36795e+01 

0.13512e+01 

0.31989e-08 

0.46930e+00 

0.12621e+02 

0.15646e+01 

0.33163e-08 

0.48941e+00 

0.26041e+03 

0.12445e+02 

0.35160e-08 

0.51098e+00 

0.23979e+03 

0.1 1604e+02 

0.36703e-08 

0.53089e+00 

0.12372e+02 

0.15638e+01 

0.37506e-08 

0.55509e+00 

0.67293e+00 

0.36566e+00 

0.37596e-08 

0.56934e+00 

0.12442e+00 

0.31993e-01 

0.37163e-08 

0.59072e+00 

0.79175e-02 

0.16584e-02 

0.35807e-08 

0.93846e+00 

0.80532e-02 

0.13858e-02 

0.62581e-08 

0.98464e+00 

0.64324e-02 

0.10296e-02 

0.70564e-08 

0.10000e+01 

0.57903e-02 

0.92241e-03 

0.41018e-04 

As  was  the  case  for  the  simpler  problem,  the  work  done  by  the  code  is  drastically  effected  by 
whether  or  not  u  is  included  in  the  estimate.  Table  14. 3j  summarizes  the  work  performed  by  the 
codes  at  various  tolerances. 


Table  14. 3j.  Work  performed  for  various  tolerances. 


TOL 

NDIFF 

No.  Steps 

0.10000e+00 

2 

61 

1 

23 

0.50000e-01 

2 

86 

1 

28 

0.25000e-01 

2 

147 

1 

34 

0.12500e-01 

2 

265 

1 

46 

0.62500e-02 

2 

289 

1 

45 

0.31250e-02 

2 

425 

1 

47 

0.15625e-02 

2 

1328 

1 

55 

0.78125e-03 

2 

1964 

1 

61 

0.39063e-03 

2 

1968 

1 

63 

Figure  14.2  shows  the  relative  improvement  in  efficiency  of  the  altered  code.  The  axes  are 
labelled  as  before,  and  results  are  perhaps  even  more  striking. 


Figure  14.2.  Accuracy  in  u  vs.  Work  Performed,  steep  gradient. 
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Problem  2 


x'(t)  =  u(t) 

u'(t )  =  f(x,u,t)  +  Gt\ 

C{x)  =0 

This  is  an  index  three  problem.  We  transform  it  to  the  following  equivalent  index  two  problem 
as  described  by  Gear  et  al.  (1984). 

x'(t)  =  u(t)  4-  Gt  ii 
u'(t)  =  f(x,u,t)  +  GtX 
C(x)  =0 
Gu  —  0 

We  assume  that  suitable  initial  values  have  been  provided.  For  example,  we  may  have  x0  and  u0 
given  such  that  the  values  satisfy  the  constraints  G(x0)  =  0  and  Gu0  —  0.  In  addition  /./0  =  0.  If 
the  other  initial  values  X0  and  u0'  are  not  known  they  will  need  to  be  computed. 

The  DAEs  lead  to  the  following  equations  for  the  errors 


< 

en 

C10  -h  0  G 

a  hp-  -h Gt  0 

dx  Ox  Ou 

—  0  0  0 

-l 

h.  T  -f-  h  if  j 
hr  +  h  //j 

'/2 

e : 

Ox 

—  G  0  0 

Ox 

<h 

We  solved  a  particular  example  of  the  above  problem.  In  this  example  we  consider  a  pendu¬ 
lum  of  unit  length  and  mass  swinging  under  the  force  of  unit  gravity.  By  using  the  constraint 
x2  -\-  y2  —  \  and  its  derivative  (xu  +  yv  =  0),  the  index  is  reduced  to  two. 

The  problem  may  be  written  as 


xu  4-  yv  =  0 


We  again  solved  the  problem  using  the  single  precision  version  of  the  code  DASSL  and  the 
results  are  summarized  here.  The  errors  in  the  computed  solution  for  various  tolerances  at  t=2n 
are  given  in  table  14.3a.  Note  that  in  this  case,  the  accuracy  of  the  altered  integrator  is  generally 
as  good  as  that  of  the  original,  for  given  tolerance. 


Table  14.4a.  Errors  in  computed  solutions  to  pendulum  problem. 


TOL 

NDIFF 

efx) 

efy) 

e(u) 

e(v) 

efX) 

max.  err. 

0.1 

6 

0.4166 

0.3271 

0.2555 

0.7126 

0.6214 

0.7126 

4 

0.3388 

0.2778 

0.1751 

0.6285 

0.4514 

0.6285 

0.01 

6 

0.0444 

0.0559 

0.0313 

0.1571 

0.1863 

0.1863 

4 

0.0404 

0.0520 

0.0430 

0.0600 

0.0039 

0.0600 

0.0001 

6 

0.0010 

0.0013 

0.0019 

0.0005 

0.0028 

0.0028 

4 

0.00007 

0.00011 

0.00058 

0.00102 

0.00993 

0.00993 

As  before,  the  amount  of  work  performed  depends  on  the  components  used  in  the  error 
estimator.  Table  14.3b  shows  the  great  improvement  in  efficiency  from  removing  the  con¬ 
trol. 


Table  14.4b.  Work  performed  at  various  tolerances. 


TOL 

NDIFF 

No.  Steps 

No.  Func.  Evals. 

0.1 

6 

27 

75 

4 

22 

61 

0.01 

6 

78 

225 

4 

47 

137 

0.0001 

6 

347 

1053 

4 

111 

327 
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